Spatial analysis of cod and flounder stomach content in relation to density using sdmTMB
# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 12))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/analysis/stomach_fullness_spatial_cache/html")
# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
ggplot(swe_coast_proj) + geom_sf()
# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 6),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Make default base map plot
plot_map_raster <-
ggplot(swe_coast_proj) +
geom_sf(size = 0.3) +
labs(x = "Longitude", y = "Latitude") +
theme_facet_map(base_size = 14)
# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
xy <- data.frame(ID = 1:length(x), X = x, Y = y)
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
return(as.data.frame(res))
}
theme_clean <- function() {
theme_minimal(base_family = "Barlow Semi Condensed") +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(family = "BarlowSemiCondensed-Bold"),
axis.title = element_text(family = "BarlowSemiCondensed-Medium"),
strip.text = element_text(family = "BarlowSemiCondensed-Bold",
size = rel(1), hjust = 0),
strip.background = element_rect(fill = "grey80", color = NA))
}
# These data are for total and prey specific stomach models
cod <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> time_period = col_character(),
#> quarter_fact = col_character(),
#> pred_weight_source = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
cod <- cod %>%
mutate(year = as.integer(year),
quarter = as.factor(quarter),
depth2_sc = depth - mean(depth),
saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
filter(year > 2014) %>%
filter(!quarter == 2) %>%
drop_na(predfle_density_sc, predcod_density_sc) %>%
droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#> converted 'quarter' from double to factor (0 new NA)
#> new variable 'depth2_sc' (double) with 102 unique values and 0% NA
#> new variable 'saduria_entomon_per_mass' (double) with 1,666 unique values and 0% NA
#> new variable 'tot_prey_biom_per_mass' (double) with 7,108 unique values and 0% NA
#> new variable 'depth_sc' (double) with 102 unique values and 0% NA
#> filter: removed 4,936 rows (58%), 3,593 rows remaining
#> filter: no rows removed
#> drop_na: no rows removed
fle <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/fle_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> pred_weight_source = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
fle <- fle %>%
mutate(year = as.integer(year),
quarter = as.factor(quarter),
depth2_sc = depth - mean(depth),
saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
filter(!quarter == 2) %>%
drop_na(predfle_density_sc, predcod_density_sc) %>%
droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#> converted 'quarter' from double to factor (0 new NA)
#> new variable 'depth2_sc' (double) with 80 unique values and 0% NA
#> new variable 'saduria_entomon_per_mass' (double) with 596 unique values and 0% NA
#> new variable 'tot_prey_biom_per_mass' (double) with 1,699 unique values and 0% NA
#> new variable 'depth_sc' (double) with 80 unique values and 0% NA
#> filter: no rows removed
#> drop_na: no rows removed
# Plot data (with smoother)
ggplot(fle, aes(predfle_density_sc, tot_prey_biom_per_mass)) +
geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3))
ggplot(fle, aes(predcod_density_sc, tot_prey_biom_per_mass)) +
geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3))
ggplot(cod, aes(predfle_density_sc, tot_prey_biom_per_mass)) +
geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3))
ggplot(cod, aes(predcod_density_sc, tot_prey_biom_per_mass)) +
geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3))
# And now read in pred_grid2 which has oxygen values at location and time and depth:
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_interactions/main/data/pred_grid2.csv")
# Clean
pred_grid2 <- pred_grid2 %>%
mutate(year = as.integer(year)) %>%
drop_na(depth)
# Add ices_rect
pred_grid2$ices_rect <- ices.rect2(pred_grid2$lon, pred_grid2$lat)
sdmTMB models with densities as covariates to stomach content data. First make mesh# Cod
pcod_spde <- make_mesh(cod, c("X", "Y"), n_knots = 150, type = "kmeans", seed = 42)
plot(pcod_spde)
pfle_spde <- make_mesh(fle, c("X", "Y"), n_knots = 70, type = "kmeans", seed = 42)
plot(pfle_spde)
# Model total prey biomass
# Cod
mcod <- sdmTMB(
data = cod,
formula = tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc,
time = "year", spatiotemporal = "IID", spatial = "on", mesh = pcod_spde, family = tweedie())
tidy(mcod)
#> term estimate std.error
#> 1 quarter1 -4.67446730 0.18369183
#> 2 quarter4 -4.62368659 0.15949784
#> 3 as.factor(year)2016 0.16511474 0.19693266
#> 4 as.factor(year)2017 0.30589646 0.20053012
#> 5 as.factor(year)2018 0.40552288 0.21650533
#> 6 as.factor(year)2019 0.27353933 0.22597431
#> 7 as.factor(year)2020 0.33937416 0.23827154
#> 8 predfle_density_sc -0.08755258 0.05474359
#> 9 predcod_density_sc 0.03191979 0.05142004
summary(mcod)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc +
#> Formula: predcod_density_sc
#> Time column: "year"
#> Mesh: pcod_spde
#> Data: cod
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> quarter1 -4.67 0.18
#> quarter4 -4.62 0.16
#> as.factor(year)2016 0.17 0.20
#> as.factor(year)2017 0.31 0.20
#> as.factor(year)2018 0.41 0.22
#> as.factor(year)2019 0.27 0.23
#> as.factor(year)2020 0.34 0.24
#> predfle_density_sc -0.09 0.05
#> predcod_density_sc 0.03 0.05
#>
#> Dispersion parameter: 0.53
#> Tweedie p: 1.74
#> Matern range: 15.22
#> Spatial SD: 0.12
#> Spatiotemporal SD: 0.61
#> ML criterion at convergence: -9904.608
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
cod$resids_mcod <- residuals(mcod) # randomized quantile residuals
qqnorm(cod$resids_mcod); abline(a = 0, b = 1)
# Flounder
mfle <- sdmTMB(
data = fle,
formula = tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc,
time = "year", spatiotemporal = "IID", spatial = "on", mesh = pfle_spde, family = tweedie())
tidy(mfle)
#> term estimate std.error
#> 1 quarter1 -5.054968610 0.27046794
#> 2 quarter4 -4.560112128 0.24061922
#> 3 as.factor(year)2016 -0.028909687 0.27451881
#> 4 as.factor(year)2017 0.301285213 0.26933035
#> 5 as.factor(year)2018 0.672111210 0.33586941
#> 6 as.factor(year)2019 0.620685733 0.28679830
#> 7 as.factor(year)2020 1.026916621 0.28115647
#> 8 predfle_density_sc 0.002925394 0.09706364
#> 9 predcod_density_sc 0.055584176 0.11765389
summary(mfle)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc +
#> Formula: predcod_density_sc
#> Time column: "year"
#> Mesh: pfle_spde
#> Data: fle
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> quarter1 -5.05 0.27
#> quarter4 -4.56 0.24
#> as.factor(year)2016 -0.03 0.27
#> as.factor(year)2017 0.30 0.27
#> as.factor(year)2018 0.67 0.34
#> as.factor(year)2019 0.62 0.29
#> as.factor(year)2020 1.03 0.28
#> predfle_density_sc 0.00 0.10
#> predcod_density_sc 0.06 0.12
#>
#> Dispersion parameter: 0.12
#> Tweedie p: 1.50
#> Matern range: 39.88
#> Spatial SD: 0.57
#> Spatiotemporal SD: 0.54
#> ML criterion at convergence: -4667.908
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
fle$resids_mfle <- residuals(mfle) # randomized quantile residuals
qqnorm(fle$resids_mfle); abline(a = 0, b = 1)
# Model Saduria contents
# Cod
mcodsad <- sdmTMB(
data = cod,
formula = saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc,
time = "year", spatiotemporal = "IID", spatial = "on", mesh = pcod_spde, family = tweedie())
tidy(mcodsad)
#> term estimate std.error
#> 1 quarter1 -8.87619104 0.7663188
#> 2 quarter4 -9.08924932 0.7214733
#> 3 as.factor(year)2016 -0.82029185 0.7820977
#> 4 as.factor(year)2017 -1.36653144 0.7993448
#> 5 as.factor(year)2018 -0.26043020 0.8378119
#> 6 as.factor(year)2019 -0.88611599 0.9053763
#> 7 as.factor(year)2020 -0.01723844 0.8681083
#> 8 predfle_density_sc -0.60227548 0.3012112
#> 9 predcod_density_sc 0.39255981 0.2620735
summary(mcodsad)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc +
#> Formula: predcod_density_sc
#> Time column: "year"
#> Mesh: pcod_spde
#> Data: cod
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> quarter1 -8.88 0.77
#> quarter4 -9.09 0.72
#> as.factor(year)2016 -0.82 0.78
#> as.factor(year)2017 -1.37 0.80
#> as.factor(year)2018 -0.26 0.84
#> as.factor(year)2019 -0.89 0.91
#> as.factor(year)2020 -0.02 0.87
#> predfle_density_sc -0.60 0.30
#> predcod_density_sc 0.39 0.26
#>
#> Dispersion parameter: 0.17
#> Tweedie p: 1.49
#> Matern range: 36.25
#> Spatial SD: 1.88
#> Spatiotemporal SD: 1.56
#> ML criterion at convergence: -1024.766
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
cod$resids_mcodsad <- residuals(mcodsad) # randomized quantile residuals
qqnorm(cod$resids_mcodsad); abline(a = 0, b = 1)
# Flounder
mflesad <- sdmTMB(
data = fle,
formula = saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc,
time = "year", spatiotemporal = "IID", spatial = "on", mesh = pfle_spde, family = tweedie())
tidy(mflesad)
#> term estimate std.error
#> 1 quarter1 -8.7008421 1.0345886
#> 2 quarter4 -7.9176420 0.9779841
#> 3 as.factor(year)2016 0.7217416 0.7990698
#> 4 as.factor(year)2017 0.5084692 0.7999960
#> 5 as.factor(year)2018 1.8324164 0.9723237
#> 6 as.factor(year)2019 0.1362148 0.8635915
#> 7 as.factor(year)2020 1.7146474 0.8095086
#> 8 predfle_density_sc -0.1720328 0.3376522
#> 9 predcod_density_sc 0.6932675 0.4896959
summary(mflesad)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc +
#> Formula: predcod_density_sc
#> Time column: "year"
#> Mesh: pfle_spde
#> Data: fle
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> quarter1 -8.70 1.03
#> quarter4 -7.92 0.98
#> as.factor(year)2016 0.72 0.80
#> as.factor(year)2017 0.51 0.80
#> as.factor(year)2018 1.83 0.97
#> as.factor(year)2019 0.14 0.86
#> as.factor(year)2020 1.71 0.81
#> predfle_density_sc -0.17 0.34
#> predcod_density_sc 0.69 0.49
#>
#> Dispersion parameter: 0.14
#> Tweedie p: 1.48
#> Matern range: 66.65
#> Spatial SD: 2.44
#> Spatiotemporal SD: 1.18
#> ML criterion at convergence: -1440.375
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
fle$resids_mflesad <- residuals(mflesad) # randomized quantile residuals
qqnorm(fle$resids_mflesad); abline(a = 0, b = 1)
mcod_ef <- tidy(mcod, effects = "fixed", conf.int = TRUE) %>%
filter(term %in% c("predfle_density_sc", "predcod_density_sc")) %>%
mutate(species = "Cod", response = "Total prey")
#> filter: removed 7 rows (78%), 2 rows remaining
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
mfle_ef <- tidy(mfle, effects = "fixed", conf.int = TRUE) %>%
filter(term %in% c("predfle_density_sc", "predcod_density_sc")) %>%
mutate(species = "Flounder", response = "Total prey")
#> filter: removed 7 rows (78%), 2 rows remaining
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
mcodsad_ef <- tidy(mcodsad, effects = "fixed", conf.int = TRUE) %>%
filter(term %in% c("predfle_density_sc", "predcod_density_sc")) %>%
mutate(species = "Cod", response = "Saduria")
#> filter: removed 7 rows (78%), 2 rows remaining
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
mflesad_ef <- tidy(mflesad, effects = "fixed", conf.int = TRUE) %>%
filter(term %in% c("predfle_density_sc", "predcod_density_sc")) %>%
mutate(species = "Flounder", response = "Saduria")
#> filter: removed 7 rows (78%), 2 rows remaining
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
plot_df <- bind_rows(mcod_ef, mfle_ef, mcodsad_ef, mflesad_ef) %>%
mutate(term = ifelse(term == "predcod_density_sc", "Cod density", "Flounder density"))
#> mutate: changed 8 values (100%) of 'term' (0 new NA)
# Plot effects
plot_df %>%
ggplot(., aes(term, estimate, color = factor(species))) +
geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.75) +
geom_point(size = 3, position = position_dodge(width = 0.2)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
position = position_dodge(width = 0.2), size = 1) +
scale_color_brewer(palette = "Dark2", name = "Species") +
labs(x = "Predictor", y = "Standardized coefficient") +
theme_light(base_size = 12) +
facet_wrap(~response) +
theme(legend.position = "bottom") +
NULL
ggsave("figures/stomach_content_effect_size.png", width = 6.5, height = 6.5, dpi = 600)
# Flounder density
nd_cod_fle <- data.frame(predfle_density_sc = seq(min(cod$predfle_density_sc),
max(cod$predfle_density_sc), length.out = 100))
nd_cod_fle$year <- 2018L
nd_cod_fle$predcod_density_sc <- 0
nd_cod_fle$depth_sc <- 0
nd_cod_fle$quarter <- factor(4)
sad_pred_cod_fle <- predict(mcodsad, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)
tot_pred_cod_fle <- predict(mcod, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)
# Cod density
nd_cod_cod <- data.frame(predcod_density_sc = seq(min(cod$predcod_density_sc),
max(cod$predcod_density_sc), length.out = 100))
nd_cod_cod$year <- 2018L
nd_cod_cod$predfle_density_sc <- 0
nd_cod_cod$depth_sc <- 0
nd_cod_cod$quarter <- factor(4)
sad_pred_cod_cod <- predict(mcodsad, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)
tot_pred_cod_cod <- predict(mcod, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)
# Flounder density
nd_fle_fle <- data.frame(predfle_density_sc = seq(min(fle$predfle_density_sc),
max(fle$predfle_density_sc), length.out = 100))
nd_fle_fle$year <- 2018L
nd_fle_fle$predcod_density_sc <- 0
nd_fle_fle$depth_sc <- 0
nd_fle_fle$quarter <- factor(4)
sad_pred_fle_fle <- predict(mflesad, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)
tot_pred_fle_fle <- predict(mfle, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)
# Cod density
nd_fle_cod <- data.frame(predcod_density_sc = seq(min(fle$predcod_density_sc),
max(fle$predcod_density_sc), length.out = 100))
nd_fle_cod$year <- 2018L
nd_fle_cod$predfle_density_sc <- 0
nd_fle_cod$depth_sc <- 0
nd_fle_cod$quarter <- factor(4)
sad_pred_fle_cod <- predict(mflesad, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)
tot_pred_fle_cod <- predict(mfle, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)
# Saduria models
p1 <- ggplot(sad_pred_fle_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "Saduria in flounder stomach [g/g]", x = "")
p2 <- ggplot(sad_pred_fle_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "", x = "")
p3 <- ggplot(sad_pred_cod_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "Saduria in cod stomach [g/g]", x = "Scaled flounder density")
p4 <- ggplot(sad_pred_cod_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "", x = "Scaled cod density")
p1 + p2 + p3 + p4 + plot_layout(ncol = 2)
ggsave("figures/marginal_effects_saduria.png", width = 6.5, height = 6.5, dpi = 600)
# Total stomach content models
p5 <- ggplot(tot_pred_fle_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "Prey biomass in flounder stomach [g/g]", x = "")
p6 <- ggplot(tot_pred_fle_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "", x = "")
p7 <- ggplot(tot_pred_cod_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "Prey biomass in cod stomach [g/g]", x = "Scaled flounder density")
p8 <- ggplot(tot_pred_cod_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "", x = "Scaled cod density")
p5 + p6 + p7 + p8 + plot_layout(ncol = 2)
ggsave("figures/marginal_effects_total.png", width = 6.5, height = 6.5, dpi = 600)
knitr::knit_exit()